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Using a numerical simulation of the classical dynamics of the plane- wave and flat 
space matrix models of M-theory, we study the thermalization, equilibrium thermo- 
dynamics and fluctuations of these models as we vary the temperature and the size 
of the matrices, N. We present our numerical implementation in detail and sev- 
eral checks of its precision and consistency. We show evidence for thermalization by 
matching the time-averaged distributions of the matrix eigenvalues to the distribu- 
tions of the appropriate Traceless Gaussian Unitary Ensemble of random matrices. 
We study the autocorrelations and power spectra for various fluctuating observables 
and observe evidence of the expected chaotic dynamics as well as a hydrodynamic 
type limit at large N, including near-equilibrium dissipation processes. These config- 
urations are holographically dual to black holes in the dual string theory or M-theory 
and we discuss how our results could be related to the corresponding supergravity 
black hole solutions. 



I. INTRODUCTION 

Some of the greatest successes of the holographic dualities have been the computations 
of transport coefficients of strongly coupled quantum field theories, beginning with the work 
of [I]. These computations replace a problem of large N (high dimensional gauge group), 
finite temperature, strongly coupled quantum gauge theory dynamics by the much simpler 
problem of analyzing classical Green's functions in various black hole backgrounds with 
boundary conditions that reflect our understanding of response theory. 

The direct calculation of the transport coefficients in the strongly coupled quantum field 
theory is currently beyond reach. The main problem is that response theory at finite fre- 
quency requires us to solve the real time evolution of an interacting quantum system of many 
degrees of freedom. If this is done in a path integral by Monte-Carlo sampling, the dynam- 
ical situation we want to consider has a severe sign problem that is currently unsolved. In 
some instances there are workarounds (for a recent review see [2]). 

The calculation in [1] and subsequent work is, to leading order and when properly nor- 
malized [3] , completely independent of N and all the precise details of the field theory. These 
are classical gravity computations in AdS black hole backgrounds. Corrections appear when 
the gravitational theory suffers quantum corrections, i.e., when the curvature is large some- 
where in Planck units, or when the truncation to gravity breaks down, as when the black 
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hole has curvature radii of order the string scale. In the dual field theory, this corresponds 
to finite iV or weak coupling. In such cases it may pay off to approach the problem from the 
other side of the duality, starting with a field theory Hamiltonian and studying its dynamics 
and then examining its relation to the gravity limit. This paper falls into this category of 
work. 

Given our inability to solve strongly coupled quantum theories directly, we could instead 
try to drop the quantum adjective and solve "strongly coupled classical theories." At first, 
this might seem impossible, but we have to consider that strong coupling in quantum theories 
also refers to the strength of the terms in the action that make the system non-linear, relative 
to the terms corresponding to a linear system Q Thus, the regime of interest in classical 
physics is that of a dynamical system in the very non-linear regime. We also want to study 
it at large N. Here we expect that the strong non-linearity will cause chaotic dynamics and 
that large iV should be understood as a thermodynamic limit, where the number of degrees 
of freedom grows as iV 2 but some aspects of the dynamics are iV-independent. Beyond 
the existence of good thermodynamic state variables, we can look for collective modes to 
emerge, akin to hydrodynamic variables, that indicate collective time dependent dynamics 
also roughly independent of N. 

Classical non-linear field theories have an infinite number of degrees of freedom and suffer 
from the ultraviolet (UV) catastrophe. The UV catastrophe is cured by reintroducing the 
Planck constant. Indeed, this is how Planck introduced the eponymous constant in the first 
place, giving birth to quantum mechanics. This would seem to stop this idea of studying 
classical non-linear field theory dynamics in its tracks. However, thinking of a field theory 
expanded in Fourier modes, a finite h freezes the dynamics of most of the modes to be in 
their ground state (or their adiabatic ground state given a configuration of the low frequency 
modes). These frozen modes are those above a cutoff, which is determined by the dynamics 
and the initial conditions. So, one is only dealing with finitely many active degrees of 
freedom and the initial problem of an infinite number of degrees of freedom in field theory 
can be solved, if we only knew the precise details of how most of the degrees of freedom 
"freeze out." 

This problem is solved without any work if we study dynamical systems with finitely 
many degrees of freedom in the first place. In that case, we don't need to address the UV 
catastrophe at all. In those systems, h — » means that we are studying the system at 
large quantum numbers (very large energies compared to the gap in systems with a discrete 
spectrum). This is a regime where typical states are well described by classical statistical 
mechanics. 

With this in mind, our main purpose is to examine the real time classical statistical me- 
chanics of certain matrix models that have finitely many degrees of freedom. This is exactly 
the same methodology as in molecular dynamics simulations (see, e.g., jl] for a review). We 
study both equilibrium configurations, with their associated equations of state, and also sim- 
ple transport processes, or more precisely, out of equilibrium relaxation. We study the latter 
via fluctuations of the appropriate variables and by invoking the fluctuation-dissipation the- 
orem. We discuss whether such classical models can be used to study holographic dualities, 
where we also have some gravitational information. That is, to what extent can the classical 
dynamics of large N matrix models encode gravitational information? More broadly, are the 

1 In perturbation theory in quantum mechanics this can also be the regime where energy denominators 
are large, but off-diagonal terms in the perturbation are even larger, so the standard conditions for 
perturbation theory to make sense do not necessarily apply. 
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dynamics compatible with our expectations from string theory, including beyond the gen- 
eral relativity/supergravity regime? We are not able to completely answer these questions, 
but do make some relevant qualitative comparisons and, especially in Sec. VII , offer some 
informed speculations on the full story, which is a subject of continued research. 

We work with the BFSS jS] and BMN [5] matrix models, which have well known holo- 
graphic dual descriptions. We simulate the real time classical dynamics of these models. In 
our simulations the generic late time behavior for the initial conditions we study are finite 
temperature, equilibrium configurations. For the reasons we discussed above, it is reason- 
able to treat these holographic matrix models classically at high temperatures, which is the 
regime of very large quantum numbers. The associated black holes have curvatures that are 
large in string units [7]. Thus, from the gravitational perspective we would be describing 
very stringy back holes with our analysis. 

We also find that there are collective dynamical variables analogous to hydrodynamic 
degrees of freedom. That is, there are certain aspects of the near equilibrium dynamics 
that are independent of the number of microscopic degrees of freedom. This was not obvi- 
ously guaranteed. Hydrodynamics usually requires a geometric coarse graining of degrees 
of freedom confined to small regions. Holographic systems are usually made of D-branes. 
The holographic degrees of freedom are the strings stretching between the branes, which are 
extended objects that generally do not localize to small volumes. In the case of AdS bulk 
geometries, even small regions on the AdS boundary have infinitely large volumes in AdS, 
so this issue does not apply directly. One does hydrodynamics on the boundary of AdS and 
not in the bulk, as in pQ . In the case of the BMN matrix model, the conformal boundary of 
the dual plane wave geometry has no spatial extent [HI E] and so one cannot have transport 
there. Hence, hydrodynamic behavior is not an obvious possibility in this case. 

However, the membrane paradigm for black holes suggests that there should be transport 
phenomena for the degrees of freedom in the near horizon limit of our black hole type objects. 
We still have to fret about the possibility of a phase transition: that the classical regime of 
the matrix model (at large N) has nothing to do with the black hole dynamics we would 
like to study. In [TO] the authors argued that the transition from branes or strings to black 
holes is a smooth process, which also suggests some form of collective behavior in the near 
horizon dynamics, but it might be very distorted from the classical gravity regime. See also 
[TTj . which also argues against any such phase transition at any finite N or finite 't Hooft 
coupling. 

We are also motivated to study this problem by the fast scrambling conjecture of [12] 
and the hope that we can find some numerical handle to study it that does not involve 
gravitational arguments in the first place. We are continuing the study that was initiated in 
[13], where some notion of fast classical scrambling was shown for the BMN matrix model. 
Another analysis in the BFSS matrix model was carried out in [TJ], which also showed 
fast classical thermalization in the BFSS matrix model for some set of initial conditions. 
However, these studies are not in the same spirit as the Sekino-Susskind setup, which seems 
to fundamentally require h ^ in its formulation. Unfortunately, we have nothing new to 
say in regards to this conjecture, even though the principal example that was argued to 
thermalize fast was exactly the BFSS matrix model. The crucial property of the model, 
shared by the ones we study, was that all of its degrees of freedom are coupled to each other 
in the Lagrangian. Indeed, it has been argued in [15] that non-commutative field theory, 
which makes the degrees of freedom more non-local, scrambles faster than the corresponding 
commutative version of the field theory dynamics, but this was again done via a gravitational 
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computation. 

Other attempts to study the fast scrambling conjecture directly include [IS] , which studies 
certain toy quantum mechanical models, but not holographic systems. See the series of 
papers [l~Tl - fT9] for a more detailed treatment of holographic fast scramblers, though these 
do not study matrix models directly. The papers [TTJ [201 EI] slightly predate the conjecture 
and study matrix model dynamics and thermodynamics perturbatively at weak-coupling, as 
well as the breakdown of perturbation theory, and so are complementary to this work. 

The paper is organized as follows. In Sec. [IT] we discuss the regime where our calcula- 
tions are valid, and we discuss the classes of observables that can be considered as well as 
a discussion on what it means for such matrix model to behave hydrodynamically. In Sec. 



Ill we describe the algorithm we use to evolve configurations. Of particular importance is 



that the dynamics requires a Gauss' law constraint, and our checks that the constraint is 
preserved to high accuracy. In Sec. IV we show that the late time behavior of the matrix 
model dynamics seems to thermalize and we present tests of this property, we also solve a 
puzzle on how to compute the temperature, where two independent measurements at first 
glance seem to disagree. In Sec. [V] we compute the power spectra of interesting observ- 
ables, especially in the BFSS matrix model. We show that the autocorrelation functions of 
interesting observables have smooth power spectrum (indicative of chaos) and have a well 
defined large N limit suggesting hydrodynamics behavior. In Sec. VI we study the large N 
factorization of the statistical quantities of interest. We show numerically that in the clas- 
sical dynamics various quantities behave to leading order as gaussians (free fields) and that 
a 1/N expansion is applicable to study various normalized correlation functions. Finally we 
conclude. 



II. OBSERVABLES AND SYMMETRY 

The objective of this paper is to study the classical dynamics of both the BFSS and BMN 
matrix models and to see to what extent we can extract lessons about holography from this 
study. Before we formulate our approach to that problem, we discuss the regime where our 
calculations are valid. 

The regime where the four-dimensional M = 4 SYM theory is dual to a semi-classical 
(super) gravity theory is large N and strong 't Hooft coupling [22]. The strong coupling 
regime implies that g\ M Nh ^ 1, so it involves h in a crucial way. One should be careful 
in interpreting this equation. In dimensions other than four, like our + 1 dimensional 
matrix models, the Yang-Mills coupling constant has units, so the left hand side can not be 
compared to the right hand side without choosing a state and multiplying by appropriate 
quantum numbers to obtain a dimensionless ratio. 

If we choose a thermal state at temperature T for an oscillator degree of freedom with 
angular frequency u>, we can ask if thermal fluctuations are larger than quantum fluctuations 
for that degree of freedom. This happens when k-^T ^> ftw. So, one can be in the small 
H regime if the temperature is high enough. For relativistic quantum field theories there is 
always some u where k-^T < hu>. Such degrees of freedom would be responsible for the UV 
catastrophe. On the other hand, for oscillators with small u the left hand side is much larger 
than the right hand side and the corresponding oscillators are at high occupation quantum 
numbers. The dynamics of these low frequency modes is controlled by classical physics. The 
classical world meets the quantum world in the intermediate regime. Roughly speaking, 
the correspondence principle in quantum mechanics should let us interpolate between the 
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classical and the quantum regimes. In the BMN and the BFSS matrix model we only have 
a finite number of degrees of freedom, so the UV catastrophe issue is avoided, but we can 
still try to push the system to the correspondence limit, in a manner we describe below. 

Typical quantum states are superpositions of position eigenstates, so if we are to match 
various physical quantities of the quantum system we should either average over positions 
or smear the classical states to a volume of h for each canonical pair of variables. If the 
system is chaotic, most energy eigenstates behave as if they are thermal for sufficiently small 
subsystems [23J (one has to make allowances if there are conserved quantities which don't 
thermalize), so the correspondence principle suggests that we should study the statistical 
properties of the thermal ensemble to study the coarse grained properties of the quantum 
states. In this paper we do not calculate anything in the regime where quantum effects start 
making a difference, but we keep in mind that in the end we want to understand the system 
in the quantum regime. 

Now, let us go back to the BMN and BFSS matrix models. The regime of interest for us 
is the large N regime, as dictated by holography. In this regime, the number of degrees of 
freedom grows like N 2 . We will see that the large N limit is not only a thermodynamic limit 
but that we also observe a kind of hydrodynamics. Here we don't mean thermodynamic limit 
in a formal, rigorous sense, but simply that we find various state variable that remain finite 
as N grows. Ideally, we should be able to show local equilibrium and transport to claim 
hydrodynamic behavior. Unfortunately, we don't know how to make such a formulation 
from first principles, as the degrees of freedom in these matrix models are essentially non- 
local. If we take one of the matrices of the matrix model, we can interpret the eigenvalues as 
positions of D-branes |24j . while the off-diagonal elements are strings stretching between the 
branes. In the classical regime we are studying, all the off-diagonal modes are excited, so it 
is hard to define local quantities that could play the role of, e.g., densities of D-particles. 

What we can do instead is add a faraway probe and determine what it sees. In the 
BFSS matrix model such a formulation leads to an effective potential for the probe. The 
effective potential can be computed from traces of the configuration [5] convolved with 
some Green's functions that decay polynomially in the distance [23]. This is a quantum 
computation where one integrates out the off-diagonal degrees of freedom connecting the 
probe to the configuration under study. This gives rise to gravitational interactions between 
general D-brane objects and gravitons (the interactions between gravitons by integrating 
off-diagonal models is part of the original formulation of the BFSS matrix model [5] 0). 
The natural candidates for hydrodynamic variables are these traces of various products of 
matrices appearing in the effective potential. They act as moments of the distributions of 
matter in the effective potential for a faraway probe. 

We need to show that the dynamics of these collective modes are roughly independent 
of N to justify calling their dynamics hydrodynamic. Note that we only make claims about 
these collective degrees of freedom in a statistical sense, just as if we were considering 
the hydrodynamic variables of a system of molecules. Since we will be studying mostly 
equilibrium configurations, all we have access to is the fluctuations of these variables. Our 
results show these fluctuations have some dynamical properties independent of N. 

Specifically, we will study time dependent correlation functions of certain single trace 
observables. Consider first a single trace operator 

O lil = tv(X i ^...), (1) 

2 The simplest one loop computation was done in [2^1 [STJ , while a two loop result was obtained in [25] . 
Higher order results require information on the wave functions of the graviton states one is scattering. 
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where [i] is a mult i- index. In the brane picture this will be a source for some gravity field 
(or more generally a closed string field). Indeed, we usually find that 



where p a would be the local source of the field (if such a notion makes sense), with its cor- 
responding spin labels. This is convolved with some polynomial function of the x, which we 
call /, which can also carry angular momentum labels. Together these would get combined 
into a multipole expansion labeled by the multi-index [i]. We can decompose the product 
into spherical harmonics, and then symmetry considerations will tell us that if the configu- 
ration is spherically symmetric the averages of objects whose spin is non-zero vanishes. For 
many interesting observables we will have that the time average should vanish 



although it does not do so configuration by configuration, but only as a time average. 
Now, if we give two such observables Oi(t), we can consider averages of the form 



where we average over a trajectory (or a collection of such trajectories with the same energy). 
The correlation function Sij(a) will describe the statistical properties of the time dependent 
correlations between the observables. Indeed, such a function will encode the fluctuation- 
dissipation information. Such correlation functions can be different from zero, even if the 
individual expectation values of the 0, vanish. This is similar to studying sound modes for 
gas in a cavity. If the individual harmonics are not excited, then their average is zero, but 
there will be thermal fluctuations. These fluctuations, when properly normalized, will have 
a good thermodynamic limit, but away from this limit there can be finite size effects that 
are sensitive to the number of particles. 

We will say the system behaves hydrodynamically if a collection of the S'y(a) properly 
normalized has a large N limit where the collection Sij(a) converges to a single function of 
a for fixed ij, and for a reasonable interval of time that is short compared to the Poincare 
recursion time, but that can be much longer than any thermalization time (or scrambling 
or relaxation time) for near equilibrium dynamics. 

This may be too narrow a definition. Consider a toy model of gas in a box with a 
somewhat random shape that is temperature dependent (like a rubber balloon filled with 
an ideal gas). We would say that the hydrodynamic behavior there is independent of the 
box. However, let us imagine that we want to study hydrodynamics by looking at the 
normal modes of sound in the box, or other such decomposition into normal modes. If 
we change the temperature, we would change the shape of the box somewhat: the added 
pressure would deform the walls of the container. This would deform the harmonics of the 
box, and the collection of harmonics of the box would be temperature dependent. Such 
changes cannot be done while preserving the spectrum (even after rescaling time). Also, if 
we change the number of particles inside the box in such a way that the pressure stays the 
same, the shape of the box would not change. However, the temperature of the gas would 
change depending on the number of particles. In such a case, the modes of sound on the 
box would be independent of the temperature only after a rescaling of time, and the ratios 
of the different frequencies would be invariant, but not the frequencies themselves. 




(2) 



(3) 



S ij (a) = {O i (t)O j (t + a)) 



(4) 
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Our systems are somewhat analogous to this. After thermalization, the matrices will 
relax to an approximately spherically symmetric configuration about the center of mass. 
This is in the absence of angular momentum for the initial conditions (we will not consider 
such initial conditions in this paper). These spherical configurations grow in size if we 
increase the temperature. In the BMN system the geometry of the plane wave in which 
the configuration is embedded acts like a box, similar to how AdS acts like a mirror. If we 
increase the temperature the configuration grows in size and the external pressure changes. 
There is also an internal pressure that makes the system want to collapse: the excitation of 
the off-diagonal modes between the branes acts like a glue that makes the system shrink. If 
these are treated as harmonic oscillators, one would expect that each such harmonic oscillator 
has an energy of k^T and that the energy stored in these configurations is independent of 
the position. However, as we move a D-particle far away from the system, the effective 
frequency of these modes goes up, and there is a corresponding shrinkage of the available 
phase space for these modes. Thus, there is an entropy cost to move a D-particle away from 
the configuration and the internal pressure to hold the system together is an entropic force. 
It has been argued that this type of entropic effect leads to the gravitational force near the 
horizon of a black hole |29|. 



On the other hand, thermal pressure makes the system expand. These two forces can 
reach an equilibrium. In the very high temperature limit we expect that the internal pressure 
dominates over the external pressure, so that the shape of the container matters less, but 
we will not be able to guarantee that the system is hydrodynamic without fine tuning: we 
would need to be able to match the box shapes between different values of N . Doing this 
carefully requires a fairly detailed understanding of the phase diagram of the system. 

All of this is much simpler to study in the BFSS matrix model. One can show that 
the classical dynamics of the BFSS model has a scaling property: a configuration at a given 
energy can be rescaled to any other energy by a rescaling of time and the matrix components. 
Thus, it does not matter at what energy we study the system as the dynamics is essentially 
the same. So we have gotten rid of the temperature-dependent shape parameters. Because 
of this, we can argue that the dynamics of the thermal system only depends on N and we 
can explore how the dynamics depends on N and no other variables. For of this reason, we 
analyze the large N limit primarily in the BFSS matrix model. 



The numerical implementation has been discussed previously in [13]. We work with a 
leapfrog algorithm and we indicate how to implement the constraints in the initial conditions. 
Here we reiterate the algorithm. 

The bosonic degrees of freedom of the BMN and BFSS matrix model are the hermitian 
matrices X l=0,1,2 and y a=1 >--> 6 and their canonical conjugates Pi and Q a . The bosonic part 
of the Hamiltonian is 




III. NUMERICAL IMPLEMENTATION 





When a ^ we have the BMN matrix model. For convenience we set a = 1 in the BMN 
case. For a = we get the BFSS matrix model. 
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We have rescaled the variables so that the classical equations of motion are independent of 
h and all the quantum mechanics is hidden in the initial conditions. We have also normalized 
the mass of X to one, i.e., we measure time by the oscillation period of one of the X modes. 

Because of the U(iV) gauge symmetry we must enforce the Gauss' law constraint: 

C = [X\P i } + [Y a ,Q a ]=0. 

To solve the equations of motion we use a leapfrog algorithm. This has the virtue of pre- 
serving the constraints. The discretized matrix equations of motion read 

dV 

X t+5t = X t + P t+ st6t , P t+ st = P t _H - — St (6) 

and similarly for the Y modes. Since we have the X 1 , Pi evaluated at different times, we 
need to be a little careful with the constraint. We define 

C{t) = \X\t), P t (t + St/2)] + \Y a (t), Q a (t + t/2)] (7) 

and will set it to zero in the initial conditions. We also define the constraint at half intervals 
to be given by 

C(t + St/2) = [X\t + St), Pi(t + St/2)} + [Y a {t + St), Q a {t + St/2)}. (8) 

To show that the constraints are satisfied notice that when we evolve the constraint by using 
the equation of motion (|6|, after one half step in t we get that 

C(t + St/2) - C(t) = Yy&X\t), Pi{t + St/2)} + ... 

i 

= J2[ p i(t + &/2), Pi(t + St/2)}St + • ■ ■ = 0, (9) 

i 

which vanishes term by term. For the second half step we need to work harder, but so long 
as V is a sum of traces of products of X and Y matrices (or functions of such traces), one 
can prove that the contribution from each such trace vanishes by summing cyclically over 
the letters making the word in the trace. Hence, C(t + St) —C(t) =0 and this tells us that 
C(t) is a constant of motion of the discrete evolution. Incidentally, the same arguments 
work for angular momentum conservation laws. Our initial conditions are those for (near) 
zero angular momentum. The only place where constraint violations might appear is from 
rounding errors, so we need to check that we don't suffer from this problem. To improve 
numerical stability we use double precision numbers. To check for numerical errors, we 
record the absolute value of the constraint CV = \ tr(C 2 )| as a check for the code. We find 
that the constraint is well satisfied for the runs we perform, so we do not need to implement 
constraint damping. The equations of motion evolve hermitian matrices into hermitian 
matrices. Truncation errors in matrix multiplication can also take us away from that locus. 
We found that we needed to enforce hermiticity of the matrices every few steps, by taking 
X — > (X + X')/2 and similarly for all other matrices. We do this every time we write the 
matrix configurations. 

The main sources of difficulty in the setup are the initial conditions. For this paper, we 
have used the following initial classical configurations: 

_ (L° n 0\ ! (L\ SxA 2 (LI Sx 2 
X ~ \0 oj' X ~ [Sx\ J' X ~ \S4 

Po =(I !)' Pl < 2 = = ^. ya = 5 V a - ( 10 ) 
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in the BMN matrix model. These are the same initial conditions that were used in [IB], 
where they are explained in detail. The important point is that the 5x and 8y are generated 
by Gaussian distributions controlled by a classical estimate of the quantum uncertainty of 
the modes which is parametrized by h. This is only used in the initial conditions and is 
similar to the standard practice in molecular dynamics simulations [I]. 

To study the BFSS matrix model, we first evolve in the BMN matrix model and some 
time later we set a = and study the further evolution of the system. Once we are in the 
BFSS matrix model, we have further dynamical information: the classical system enjoys a 
scale invariance. Thus, results at one energy are completely equivalent to results at higher or 
lower energies (there is no temperature dependence on quantities, other than those dictated 
by scaling) and the only variable we have is N. This makes it easy to compare different 
values of N to each other after proper rescaling. In this situation we can test convergence 
of quantities as we increase N in a temperature independent way. 

We store the full configurations of the matrices every few steps in 5t (for the data sets 
we present here we set this number to ten unless otherwise stated), and we store other 
information for faster processing at different intervals. This is especially important for long 
simulations. We will call machine time the total run in the simulation in units of the smallest 
time step that is recorded. 

Our results for the constraint violation can be seen in Fig. [I] As the constraints tech- 



CV NCV 




FIG. 1. Constraint violation as a function of machine time and normalized constraint violation. As 
can be seen the constraint violation stays very low through the computation. The graph indicates 
values of N = 4, 7, 10, 13, 18, 27, 47, and after normalization it shows larger fluctuations for smaller 
N. These are due to the statistical fluctuations of tr(X 2 ) tr(P 2 ), which are an integral part of 
the dynamics. In the plot to the right we have selected larger intervals in machine time to aid 
visualizing the different values of N. 

nically have units, we need to normalize them, and we define the normalized constraint 
violation to be given by NCV = -iVtr(C 2 )/(tr(X 2 ) tr(P 2 ))(t) (we randomly chose X 2 , Pi) 
and N being the size of the matrices. In the simulations depicted in the figure [I] we set 
a = after machine time t = 2000. 



IV. THERMALIZATION 



In [IB] it was shown that the initial conditions ( 10 ) generate configurations of eigenvalues 
which coalesce into a uniformly oscillating blob, for example see Fig. [2j The system was 
argued to have thermalized in that time averaged distributions of the momentum degrees 
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of freedom follow a Gibbs ensemble dP dQ exp(—/3H) for some inverse temperature 0. The 
Gibbs distribution factors into a product of Gaussians because both the BFSS and the BMN 
Hamiltonians are quadratic in the momenta. It was shown that the binned eigenvalues 
collected over time follow the semicircular distribution for random matrices, as is expected 
for random Gaussian matrices in the large N limit. Here we explore in more detail the 
nature of the thermalization, the appropriate finite N Gibbs distribution, and the design of 
a correctly calibrated thermometer. 

One thing to note is that the equations of motion for the traces of the coordinates and 
their momenta are those of a harmonic oscillator: the trace of the X 1 and P l oscillate with 
period 2ir while the Y a and Q a oscillate with period Air. Although this property of the BMN 
system can be used to generate a clock for the system (as in Fig. [2]), it is undesired here. The 
trace is a protected degree of freedom that does not thermalize. To describe the thermalized 
system in a statistical manner, we must remove the trace degree of freedom from our Gibbs 
distribution. The matrices are Hermitian, and so we must also enforce that on our Gibbs 
distribution. What we are left with is the Traceless Gaussian Unitary Ensemble (TGUE), 
a means by which to select traceless random Hermitian matrices. The trace of the matrices 
represents the center of mass motion of the system and thus our partition function really 
only describes the internal degrees of freedom that can thermalize. 

In order to study ensemble quantities of the system, we must coarse grain the dynamics 
leaving only gauge invariant quantities to be studied. We can study our Gibbs distribution 
in a gauge invariant manner by focusing on eignvalues and traces. Integrating over the 
unitary degrees of freedom gives the joint probability distribution for the eigenvalues. This 
result is well known for the GUE and is simple to modify for the traceless case in which 
we are interested. The trace is invariant under a unitary transformation, and thus we 
may enforce tracelessness by inserting a delta function without affecting the removal of the 




FIG. 2. Eigenvalues of X° as a function of time for a simulation of rank 8 matrices with St = 0.001, 
v = 20.0, and h = 0.001. The abscissa measures discrete time units between recordings. The 
sinusoidal curve at the bottom is the trace of X° which serves as a clock due to its equations of 
motion. 
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FIG. 3. In the left figure we have zoomed in on the time interval [370, 470] in Figure [2j In the right 
figure we have blown up the small rectangle in the left figure to observe more closely a crossing 
between the two largest eigenvalues. The sampling rate was increased 50 fold to obtain the right 
figure, however the time scales between the two figures have been kept in sync to avoid confusion. 



unitary degrees of freedom. We define the following partition function as the integral of the 
joint probability of the eigenvalues for the TGUE 



Z\= c? JV A5(tr(EjAj)) J] lAi-A/exp -^A?j (11) 

l<i<j<N \ i=l J 

where N is the rank of the matrices of interest and /3 is a parameter analogous to the 
standard deviation for normal distributions. For us (5 is physically the inverse temperature 
of our system. Note that the polynomial in the integrand is the square of the Vandermonde 
determinant of the eigenvalues. 

Before moving on to the thermodynamics, we would like to point out something about 
the dynamics of the eigenvalues. We can move the Vandermonde determinant into the 
exponential 

Z x = /^A5(tr(EAi))exp(-^A l 2 + 2^1og|A i -A,|) (12) 

The exponential describes a quadratic potential with a logarithmic term. Thus the eigenval- 
ues should always repel each other and do not cross. Although the rogue eigenvalue in Fig. 
[2] appears to pass through the others at early times, it actually just transfers energy to the 
adjacent eigenvalue, like a Newton's cradle. This behavior can be realized in plots like Fig. 
[2] with a sufficiently small time step, a large enough sampling rate, and enough zooming. 
An example of such is shown in Fig. [3] 

To determine if our system has thermalized, the first step is to match our eigenvalue 
distribution to the distribution predicted by the partition function Z\. The Vandermonde 
determinant becomes exponentially complex with increasing iV and so a direct comparison 
is time expensive. Instead we use the probability distribution obtained by integrating the 
partition function over all but one eigenvalue, i.e., the eigenvalue probability density or level 
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FIG. 4. A histogram of the eigenvalues of P° for a simulation of rank 6 matrices sampled after 
thermalization. A total of 4 x 10 5 configurations were sampled meaning we fitted using 2.4 x 10 6 
eigenvalues. By computing the second moment, we can fit the distribution to the level density of 
the TGUE. The level density has been normalized to the total number of samples times the bin 
size. The fit has an R 2 value of 0.999904. 



density function. An explicit form of the level density for the TGUE for arbitrary N has 
recently been found in |30j . Fig. [4] shows that the eigenvalues of the traceless momentum 
matrices sampled over time after thermalization do indeed fit the predicted function. 

In order to make these comparisons, we need to know the temperature 0. If we consider 
the tracelessness of the matrices as a constraint, then we can apply the proof in Appendix 
[Alto a single traceless matrix and obtain 



T 



(tr(P 2 ))o 
N 2 - 1 



(13) 



where P is any single momentum matrix and the zero subscript indicates we take the ex- 
pectation value with respect to the TGUE. For numerical measurements, the zero subscript 
indicates averaging only over the traceless part of the matrices. The momentum contribu- 
tion to the Hamiltonian is invariant under an SO(9) transformation and one would suspect 
that we could use any momentum matrix and obtain the same temperature. The numerical 
data of Table H] shows that this is not the case. 

What we have forgotten about in developing our thermometer are the constraints. Indeed, 
if we read directly from the table under our naive assumptions, we would conclude that the 
system has two different temperatures and therefore has not thermalized, even though it 
has isotropized along the relevant directions. It would seem that either the system is not 
thermal or the thermometers are broken. To answer this puzzle notice that the naive Gibbs 
ensemble is over all P, Q matrices, but our true ensemble is over P, Q, X, Y matrices, 
subject to the gauge constraint. The gauge constraint can not be ignored: it is an essential 
part of the dynamics. For example, if we ignored it we would get the wrong specific heat 
for the system. In quantum systems this is usually cured by adding ghosts, but in this case 
we need to add the constraint to the ensemble. 
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The gauge constraint induces symmetry breaking between the P, Q variables which breaks 
the SO(9) symmetry, but preserves the SO(3) x SO(6) subgroup. The SO(9) symmetry is 
broken between the X and Y matrices by the dynamics and thus the typical values of the 
X and Y matrices are different. This breaking of S0(9) therefore induces an asymmetry in 
the Gibbs ensemble for the P, Q matrices because they appear with different coefficients in 
the constraint. We also have conservation of angular momentum in the X — P and Y — Q 
planes, so we should also remove the angular momentum degree of freedom as a constraint: 
we chose initial configurations with zero angular momentum. Thus the SO (3) x SO (6) 
subgroup is restored on averaging throughout the entire simulation and explains why in 
Table [I] it appears that each of the Pi yield the same temperature as do each of the Q a , but 
independently they seem to be different. 

Putting everything together, our true Gibbs distribution is the GUE for nine matrices 
with the constraint of tracelessness, the gauge constraint, and the total amount of angular 
momentum set to zero. Each of these constraints is linear in the momenta and thus we may 
apply the result in Appendix [A] which is a generalized equipartition theorem, to obtain an 
absolute normalization of the temperature. 

E(tr(P l 2 ))o + £<tr(Q 2 J)o = - !) - - !) - _ ^r) T = ^ " 26)T - 

i=0 o=l ^ ' 

(14) 

The lowest rank matrix which can exhibit thermalization behavior nontrivially is N = 2 and 
so we do not run into an issue of negative temperature (this would indicate that we have 
not taken into account relations between the constraints). 

We would also like to ensure that temperature measurements are independent of the 
time step parameters in our numerical implementation, that is, St. The BMN matrix model 
exhibits chaos and so shrinking the time step does not cause the numerical solution to 
converge to a solution of the equations of motion as small differences grow exponentially 
for large times. We still expect that each of the trajectories computed this way would 
lead to the same ensemble since we should be sampling the phase space according to the 
dynamically invariant measure, in a manner typical of numerical simulations of chaotic 
dynamical systems. 

In order to measure the temperature we need to measure (tr(P 2 ))o for all momenta 
matrices. We can determine the accuracy of our measurements by dividing the configurations 
into groups to obtain several sample measurements of (tr(P 2 ))o- The expectation value is 
independent of the grouping of configurations due to its linearity, but now we can obtain 



N 


(tr(P 2 ))o 


(tr(P 1 2 ))o 


(tr(P|))o 


(tr(Q?))o 


(tr(Q|))o 


<tr(Q!))o 


(tr(QD)o 


(tr(QD)o 


<tr(Q!))o 


4 
11 
23 


23.2 ±0.6 
26.9 ±0.3 
32.2 ±0.3 


23.3 ±0.4 
27.2 ±0.2 
32.2 ±0.2 


23.2 ±0.5 

27.0 ±0.3 

32.1 ±0.2 


21.3 ±0.5 
26.6 ±0.2 
31.9 ±0.2 


21.3 ±0.5 
26.5 ±0.3 
31.9 ±0.2 


21.2 ±0.6 
26.6 ±0.2 
31.9 ±0.2 


21.2 ±0.4 
26.6 ±0.3 
31.9 ±0.2 


21.3 ±0.4 
26.6 ±0.2 
31.9 ±0.2 


21.0 ±0.4 
26.5 ±0.2 
32.0 ±0.2 



TABLE I. Time average samples of the trace of the square of the traceless momenta for various 
N with v = 20.0, h = 0.001, and 5t = 0.001. For each expectation value, 20 samples were used 
but the number of configurations varies with N due to limited hard disk space. The differences 
between the values in the disjoint groups of (tr(P 2 ))o and (tr(Q 2 ))o are smaller compared to the 
differences of the values in the union of these two groups indicating the breaking of the SO (9) 
symmetry between the momenta that is present the Hamiltonian. 
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N 


St 


4 


14 


23 





005 


1 


932 ± 0.014 





16235 ± 


00050 


0.06845 ± 0.00018 





003125 


1 


932 ±0.014 





16234 ± 


00043 


0.06844 ± 0.00017 





0025 


1 


932 ±0.014 





16235 ± 


00047 


0.06843 ± 0.00016 





002 


1 


932 ± 0.016 





16234 ± 


00048 


0.06843 ± 0.00018 





00125 


1 


932 ± 0.013 





16234 ± 


00050 


0.06843 ± 0.00016 





001 


1 


932 ±0.014 





16234 ± 


00048 


0.06843 ± 0.00016 





000625 


1 


932 ± 0.015 





16234 ± 


00053 


0.06843 ± 0.00016 





0005 


1 


932 ±0.017 





16234 ± 


00051 


0.06843 ± 0.00017 





0004 


1 


932 ±0.015 





16234 ± 


00051 


0.06843 ± 0.00015 



TABLE II. Measured temperatures of thermalized system for rank 4, 14, and 23 matrices for several 
St using the same initial conditions for each N. The sampling rate is chosen such that the time 
separation between recorded configurations is kept constant, in particular (sampling rate) x (St) = 
0.05. A total of 20 samples were used for each measurement, however, the number of configurations 
per sample was decreased with increasing N. 



a standard deviation. Consecutive configurations are correlated, so some care must be 
taken when grouping configurations to make a sample. To minimize correlations between 
samples, we group configurations consecutively. Each sample will be correlated with itself, 
but different samples will only be correlated at their boundaries. This sampling process 
provides a way to measure (tr(P 2 )) with some degree of accuracy. The temperature is 



computed using equation (14) and the standard deviation is computed by summing the 
standard deviations of the (tr(P 2 )) in quadrature. 

Table [III lists the temperatures of simulations for various N and various St with the same 
initial conditions for each N . The temperatures are equal to several significant figures and 
the error bars intersect a common average value. Furthermore we find that the coefficients 
of variation are less than 1%. To claim our comparison among different St is reasonable, 
the sampling rate is chosen such that the time between recorded configurations is constant. 
We conclude that the temperature is a well defined quantity regardless of how far apart 
trajectories in phase space become due to changing the time step. 

Some simple observables that are also natural thermodynamic variables are the sizes of 
the distributions of X and Y matrix eigenvalues. We can determine how they scale with 
the temperature and N using the virial theorem. For our case, a virial can be computed for 
each X and Y matrix. Consider, for simplicity, the expression 

j t tr(X*P) = tr(Pf) + tr(X*P) = tr(P 2 ) - tr(X*d x .V), (15) 

where we are only using the traceless parts of the matrices (we subtract the trace modes, 
which are decoupled) and we don't sum over i. We then integrate over a period of time r 
and average and obtain 

^(tr(X*P)(r) - tr(X*P)(0)) = <tr(P 2 ) - tr(X i d x .V)) r (16) 

If the trajectories are bounded then the left hand side asymptotes to zero as r — > oo and we 
obtain a relation between the kinetic energy and various derivatives of the potential energy. 
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This is simplest in the BFSS matrix model. We find after summing over the X{ that 



E tr (^ 2 ) + E tr ([ x4 ' XJ ][ x4 ' XJ ]) = °- 



(17) 



4E pot . We have 



This is, we find that the kinetic energy is twice the potential energy 2E] iin 
already argued that the left hand side grows like (8iV 2 — 26)T, so we find that the right hand 
side takes the same value. The total energy in that case is 



E tot =-(8N 2 -26)T 



(18) 



At large N, we get that the energy as a function of the temperature is |(8iV 2 )T. The specific 
heat is essentially the same as that of 6N 2 harmonic oscillators and is constant. Notice that 
this result also matches the Monte-Carlo lattice simulations in the matrix model, as seen in 
Fig. 3 of |31j . Other such simulations |32j do not cover the high temperature regime. 

Another means of getting at the size of the X and Y matrices is to look at the distri- 
bution of the eigenvalues. The elements of any single coordinate matrix appear at most 
quadratically in the Hamiltonian. Integrating the Gibbs distribution exp(— (3H) over all 
momenta and all but one coordinate matrix will be a Gaussian distribution in the remaining 
coordinate matrix elements. The standard deviation will be modified due to the constraints, 
but the form of the integrand will remain unchanged. The only constraint left over is the 
tracelessness of the matrices. Thus we expect that the eigenvalues follow the level density 
of the TGUE. Observing Fig. [5] this is exactly what we see. 





2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.1 



-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 



FIG. 5. A histogram of the eigenvalues of X 2 and Y 6 for a simulation of rank 6 matrices sampled 
after thermalization. A total of 4 x 10 5 configurations were sampled meaning we fitted using 
2.4 x 10 6 eigenvalues. Each fit has an R 2 value of 0.99925 and 0.9988 respectively. 



We can also estimate the commutator squared term if we assume that the X and Y 
matrices are random. If the eigenvalues of X are of order a (we call this the size of the 
matrix), then the eigenvalues of [X 1 , X^][X l , X J ] grow like a 4 , and we get that Na 4 ~ N 2 T. 
Thus the size of the matrices grows like a ~ A^ 1 / 4 T 1//4 . This is also true for the BMN matrix 
model at high temperature. In that case, the cubic and quadratic terms in the potential 
are subleading when the size of the matrices gets large and one asymptotically matches the 
BFSS matrix model. A test of our prediction is shown in Fig. [6]with remarkable agreement. 
Our claims hold for large N and so we see larger deviations for smaller N. 
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FIG. 6. A plot of the size of X°, measured as a = (tr(X 2 )) /N, as a function of temperature for 
various N. The error bars for larger N are smaller than the point size and thus may not be visible 
in the plot. The lines are given by a = cN^-'^T 1 '^ for constant N. Doing a least squares fit gives 
the constant c = 0.504. Plots for the other coordinate matrices give identical results. 

V. POWER SPECTRA AND CLASSICAL CHAOS 

As we have seen, there is evidence for thermalization in the BMN matrix model. Similar 
considerations show that the BFSS matrix model thermalizes (this has been studied for 
different initial conditions in [H]). This should not be surprising. Both the BFSS and BMN 
matrix models result from dimensional reduction of SYM to constant configurations on either 
flat space or a sphere. It turns out that the dynamics of translation invariant configurations 
of Yang Mills theories generally exhibit chaos [331 131] and therefore the BFSS matrix model 
exhibits chaos (this was reiterated in [33]). Because of classical scale invariance of the Yang 
Mills action, this extends all the way to infinitesimal configurations of the fields. Chaos is 
also present if a mass term is added [36] , but to access the chaotic region requires finite field 
configurations. The BMN matrix model is effectively a massive version of the BFSS matrix 
model, so it should also exhibit chaos for field configurations where the fields are sufficiently 
large, but can display integrable behavior for small oscillations around vacuum states. 

In this section we analyze the chaos in both the BFSS and BMN matrix models and show 
how we can use this information to study holography. For this purpose, let us pretend that 
a configuration in our classical system represents a thermal equilibrium state in a quantum 
system. Then we would be interested in various response functions and correlation functions 
of observables in order to understand the dynamics of the thermal state. 

For example, a typical gauge invariant observable would be a trace. The simplest traces 
are those of the matrices X 1 and Pj. However, we don't gain much from studying these 
as they are decoupled and either work as a harmonic oscillator (this is the center of mass 
motion in the BMN matrix model), or they give free non-relativistic motion on flat space 
(this is the center of mass motion in the BFSS matrix model). So we need to look for 
traces of more complicated composite objects. Here it pays to notice that we are studying 
configurations with zero angular momentum. For black holes, this means we are looking 
for spherically symmetric configurations in ten dimensions, so the perturbation modes of 
the black holes may be characterized by their angular momenta. The matrices X % form a 
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9 of SO (9) in the BFSS matrix model (in the BMN matrix model the appropriate group is 
SO(6) x SO(3), which is only slightly more complicated to analyze), so it is convenient to 
study configurations that are highest weight states of SO (9) multiplets. Spherical symmetry 
then predicts that the one point correlation functions of SO (9) non-singlets in a thermal 
state would be zero, but two point functions could be non-zero if there is a singlet in the 
tensor product of the two SO (9) representations. 

Let Z = X 1 + iX 2 . This is a highest weight state of SO(9), so we can also take operators 
Ol = tr(Z L ). These will be highest weight states of symmetric traceless combinations of the 
X with angular momentum L. In the case of M = 4 SYM theory such modes constructed 
from scalars at zero temperature are protected states [37] and they are the dual description of 
gravity excitations of AdS space. These states already display incompressible and dissipation 
free hydrodynamic behavior, at zero temperature, in the form of a quantum hall droplet [38J. 
Collective excited states can be put in one to one correspondence with gravity states [39J and 
the shape of the gravity configuration is directly determined by the expectation values of 
these traces. We expect that these simplest traces are also closely related to gravity modes 
of a black hole in the in the dual of the M = 4 SYM theory at finite temperature. This leads 
us to expect that these dynamical variables are also related to gravity modes in both the 
BFSS and BMN matrix models at finite temperature. For example, they could describe how 
gravitational or dilaton partial waves (see, e.g., [IQJ HI]) are absorbed or emitted from such 
systems. In any case, these variables are important for understanding these configurations 
in detail and could help us to ultimately learn about emergent black hole phenomena like 
Hawking radiation, or whatever corresponds to it in the regime we are studying. 

The conjugate operators are Ol = ti(Z L ). An interesting two point function is then 



which, owing to time translation invariance of the thermal density matrix and the Hamil- 
tonian, can only depend on the combination (t — t'). One can also study this correlation 
in Fourier space, so that Sl(u) = J Sl(cl) exp(—ioja)da, which is how frequency dependent 
transport coefficients are usually defined. Closely related quantities can be calculated in 
gravitational setups by using the holographic dictionary in a perturbed black hole geometry 
with infalling boundary conditions at the horizon. This ultimately gives a relation between 
the quasinormal modes of asymptotically AdS black holes and CFT response functions (see 
[521454] for reviews). We need to understand how we can compute these quantities in the clas- 
sical dynamics. The key observation is that in quantum chaotic regimes we expect S^t — t') 
to be roughly independent of the microstate that we choose, even if it is a pure state, so long 
as it is a typical state of the thermal system |^J Indeed, we expect most correlation functions 
of energy eigenstates \Ei) to be approximately thermal: 



and the dependence on t — t' is guaranteed by time translation invariance of the matrix 
elements of the energy eigenstates. Here thermalization means that 5" decays to zero (usually 
exponentially), with some thermalization time r, and that the left hand side approximates 
the right hand side on time scales that are short compared to Poincare recurrence times. 



given by 




(19) 




(20) 



3 This expectation is an extension of the idea that all energy eigenstate behave as if they are thermal states 
for time independent questions |23j . Some evidence of this behavior can be found in examples |42j . 
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We extend this to more typical states, which are superpositions of energy eigenstates 
around some energy value, by averaging over t keeping t — t' fixed. We then expect 

((iP\(D L (t)6 L (t + a)\iP}) t ~ S L (a), (21) 

where is some typical state and and we average over time. The quantum theory should 
match the classical theory for these correlation functions in the correspondence limit: at 
very large quantum numbers the answer in the classical and quantum theory should be very 
similar so long as the time scales involved are relatively small compared to the Poincare 
recurrence time. For a many body system like our large N matrix models, all the time 
scales we consider will be much smaller than the recurrence time. The is an example of 
using classical statistical mechanics as an approximation to quantum statistical mechanics. 



We compute the left hand side of (21) by using a typical state of the microcanonical 
ensemble and averaging over its trajectory. We obtain these from our simulations by first 
waiting until the system thermalizes then averaging over various configurations. The func- 
tions Sl(cl) are the autocorrelation functions of the system. Let us first consider the time 
series of Oi,{t) for some L after thermalization. We display this in Fig. [7j From the figure we 
see oscillations of a typical frequency, but they are not regular nor centered on zero. Rather, 
they appear to be superposed on waves of a much longer period than the time period shown. 
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FIG. 7. Time series of the function Re(tr(Z 2 (i))). We show four different time series obtained by 
rotations of the Z by SO (9) action. They all look similar, showing approximate rotation invariance 
of the time average. We showcase the discrete data we have in one of them, so show that the time 
dependent features are well covered by our time slicing. The sample shown here is for 18 x 18 
matrices. 

To extract information from the time series, we compute Sl{oo) by taking the Fourier 
transform of O^t) and averaging over many configurations. In our case, rotation invariance 
tells us that we did nothing special by choosing (X 1 + iX 2 ) as our highest weight state. Since 
all tensor representations of SO(n) are real (X 1 — iX 2 ) can be obtained from a rotation of 



(X 1 -MX 2 ). This means that the right hand side of (21 ) only depends on the absolute value 
of a and therefore the power spectrum P{oS) = S(\u\) is an even real function of u. Thus we 
only have to display the answer for u > 0. We show power spectra for L = 2, 3, 4 in Fig. |8j 
The first observation we can make from these plots is that the power spectra are those of 
a chaotic system. If a system is integrable, we expect the system to be solvable in terms of 
action-angle variables. The angle variables are multivalued, with period 2ir or 1 depending 
on conventions, but they have very simple time dependence <fi a (t) = a (O) + u a t. Any 
single valued function on phase space can then be represented by its Fourier series in the 
angle variables and its time dependence is that of a quasi-periodic function of time, with 
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FIG. 8. The power spectrum of 5(a) = (trpf 1 + iX 2 ) L (t) trpf 1 - iX 2 ) L (t + a)) for various 
L in random units. Results shown are for 13 x 13 matrices in the BFSS matrix model after 
thermalization. The results are averaged over 15 runs of the same length, taken from splitting a 
single time series in 15 equal parts. The jiggling of the data should be interpreted as an estimate 
of the statistical error bars for each frequency. 



characteristic frequencies determined by all integer linear combinations of the u a . Thus, the 
power spectrum of its time series should display delta-function peaks at the characteristic 
frequencies of the system. The series we observe has no delta function peaks, rather it seems 
to be described by broad band noise. This is one of the standard criteria to distinguish 
chaotic from non-chaotic systems [43J. 

Now, we can also make more sense of what we see in Fig. [7j Notice that for L = 2 
there seem to be two peaks, one near zero and another one at a characteristic frequency 
ojq. Thus we may describe the signal approximately as oscillating with some characteristic 
frequency while riding on a very low frequency envelope. We also have information on other 
modes. For L = 4 we observe broader peaks located roughly in the same places as well as 
a frequency doubling of the Uo peak. For L = 3 we observe peaks at 'half period' spacings 
relative to L — 2. The reader may have noticed that in Fig. [7] we used 18 x 18 matrices, 
whereas in Fig. [8] we studied instead the system of 13 x 13 matrices and may be worried. 
The natural way to understand this is to look at how the power spectrum depends on N, 
the size of the matrices. This is shown in Fig. [9j 

We see from Fig. [9] that the power spectrum of all L = 2 modes for various values of 
N are actually very similar to each other. Each has a large peak at zero, which is more 
noticeable in a log-linear plot. Also, we see a second peak at some characteristic frequency 
which depends on the energy of the system and N. We compare various values of TV by 
finding the location of the peak and rescaling the power to make the plots lie on top of each 
other. To find the location of the peak we do a fit of the log(P(co>)) to a quadratic function 
of u in a small interval around the visual maximum. We then extract the value of u that 
corresponds to the maximum and we scale each axis of frequency to the corresponding lon 
found for each N. The main systematic error comes from the choice of the interval. We 
show the result in Fig. [10} The data have collapsed to a single graph. 



As shown suggestively in Fig. 10, the logarithm of the power spectrum seems to have 



rather distinct features that are characterized by straight lines. We can numerically compare 
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P(OJ) 




FIG. 9. The power spectrum of tr(X 1 A-iX 2 ) 2 (t) for various sizes of N x N matrices in random units 
(its Fourier transform is S(a)). Results shown are for the BFSS matrix model after thermalization. 
The results are averaged over 15 runs of the same length, taken from splitting a single time series in 
15 equal parts. We also average further over 8 rotations of the variables to increase the statistics. 
The jiggling in the curves gives a measure of the statistical error bars of the data sets. 




FIG. 10. The power spectrum of tr(X x + iX 2 ) 2 (t) for various sizes of iV~ x N matrices. The axis 
of frequency has been rescaled for each N, to the frequency ljn, and we have also rescaled the 
power spectrum. The reference frequency for each ./V is located at 1 in the graph. Results shown 
for N = 7, 10, 47. We also have drawn additional suggestive straight lines superposed on the graph 
that serve as distinctive features of the power spectrum. 



the different values of iV to get an idea of how closely the curves match by considering the 
width of the peak near zero, relative to ujv. The dimensionless width can be parametrized 
by the slope 

'd\og(P(u) N ) 



doj 



-u N 



(22) 



which is evaluated near zero and with a cutoff slightly below 0.4cuAr. Larger 7^ corresponds 
to a larger width. This is a dimensionless number that can be used to quantify how close 



the curves at different N are to each other. We show this in Fig. [TTJ As seen in the figure, 
all values of iV > 4 have close to same behavior and the differences are controlled by the 
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7n 



FIG. 11. 7n versus N. The error bars indicate systematic errors from choosing the fitting intervals. 
A fit to a single number has been done ignoring N = 4 which is an outlier by inspection. All of 
the different values of N > 4 align within the systematic errors. 

systematics of the fit, which is dominated by the choice of interval over which we compute 
the slope. This matching is necessary to have a well defined large A" limit for these time 
dependent correlation functions. We have checked that the graphs are very similar for other 
simple operators and so they all seem to have a good large A" limit. 

How should we interpret these results? One way is to conclude that the system is be- 
having hydrodynamically: there are some large A" collective variables whose time dependent 
characteristics are independent of N, up to some rescalings of the variables by the natural 
frequency of the dynamics. We just checked the simplest angular momentum mode with 
L = 2, but we can do the analysis for tr(Z L ) for various L. The plot of the power spectrum 
for various L can be seen in Fig. [12} As the reader can see, the patterns observed for low 
L in Fig. [8] persist. Notice that the logarithmic scaling of the power spectrum makes the 
pattern more regular. To investigate this in more detail we need to address further how the 




FIG. 12. Power spectrum in arbitrary units for Ol, with L = 2, . . . 10, with values of L increas- 
ing from bottom to top in the graph. The plots are vertically separated so they can be easily 
distinguished. For each L we show two such sets. These data are from N = 27. 

different A" are related to each other. We will do this in the next section. 
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We can also look at what happens when we deform the system from the BFSS matrix 
model to the BMN matrix model. It is interesting to see how symmetry breaking is imple- 



mented in the power spectra. This is depicted in Fig. 13 



P(CJ) 




FIG. 13. Power spectra in arbitrary units for O2 for various matrix combinations. The plots are 
artificially separated so they can be easily distinguished. In most of these plots the net normaliza- 
tion of tr(X 2 ) ~ tr(T 2 ) is very close to each other, as shown previously. These arise from 40 x 40 
matrices with initial velocity in our initial conditions set to v = 100. 



As shown in the figure, the power spectrum of tr(X x + iX 2 ) 2 acquires a new bump near 
where the power spectra of the BFSS matrix model had a local minimum. The bump is less 
pronounced in the tr(XY) channel and seems absent in the ti(Y l +iY 2 ) 2 channel. There is 
also a deformation of the tr(XY~) bump near the second minimum. This shows that three 
objects that had the same symmetry properties in the BFSS matrix model exhibit the broken 
symmetry of the BMN matrix model in their dynamics. If we zoom in near the bump at 
zero, we can also see small differences. The size of these bumps depends on the strength 
of the mass term in the BMN matrix model. A full analysis of such deformations would 
consider the mixing between modes in different symmetry classes and would provide some 
understanding of response theory beyond linear response: we made a finite deformation of 
the lagrangian and the response to the deformation can be measured in dynamical quantities. 

Understanding how the bump size depends on N and the effective mass, keeping the 
temperature fixed, can give us a better understanding of the phase diagram of the BMN 
matrix model and can let us fine tune the system to obtain an appropriate large N limit. 
This requires the effects of the mass term and cubic term deformations to be compatible 
with the large N scaling we obtained in Sec. IV We can analyze this in terms of their 



expected contributions to the free energy. We expect the corrections to the free energy from 
the mass terms to be of order N/j^N^-^T 1 ^ 2 , as compared with N 2 T. If we want the ratios 
of these two contributions to the free energy to stay fixed (so that we get a proper large N 
counting of the free energy), we need to scale fi 2 ~ N 1 / 2 ^ 2 for some reference temperature 
T . Performing such an analysis is beyond the scope of the present paper. 
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VI. FACTORIZATION 



A crucial aspect of large iV physics is factorization. This states that correlators have 
a large N expansion in powers of 1/N, where the leading power of iV arises from planar 
diagrams jH] and subleading corrections arise from higher genus Feynman diagrams. For 
the simplest observables, the leading expectation value of a product of observables is the 
product of the expectation values, so long as these expectation values do not vanish in 
the first place. The arguments of planarity, or more precisely, that many of the simplest 
excitations in such systems give rise to an approximately free theory, with interactions 
governed by 1/N or 1/yN corrections is an integral part of gravitational holography [22] 
(see |45j for a nice description of this physics). 

The classical dynamics we have factorizes in a trivial sense: the value of a product of 
any set of observables at time t is the product of the values. However, what we would like 
to check is that expectation values that are averaged over time have this property as well: 
that the simple degrees of freedom can be converted into approximately 'free' constituents. 
The simplest way to think about this is that the thermal fluctuations in observables are 
independent of each other for the factorized degrees of freedom, and because the degrees of 
freedom are approximately free, the fluctuations should be gaussian. In the quantum theory 
near the vacuum, there is a standard way to understand that this leads to a consistent large 
N classical dynamics [IB] . Here we want to check that there are also classical thermodynamic 
(or more precisely hydrodynamic) variables on which one can do a similar type of analysis. 

Imagine that the simulations we are doing with time evolution in the BFSS or BMN 
matrix model can be reinterpreted as a matrix model calculation 

^ f Q(X)eM- PV(X))mm 
{U[t))t - JeM-mX))uu (3) 

If the right hand side factorizes, then so does the left hand side. Indeed, the algorithm we 
are following would correspond to a hydrodynamic Monte Carlo code to compute the right 
hand side. So long as the trajectories in the system we have are sufficiently mixing, then 
the left hand side and the right hand side should match for a long enough t. 



What should be interesting to notice is that the right hand side of Eq. (23) in the 
BFSS matrix model has no quadratic term. Hence the usual argumentation based on planar 
diagrams does not hold, as there are no quadratic terms in the X variables. One can also 
argue that in the BMN matrix model at large (3 the quadratic terms matter very little, 
and it is instead the quartic term that dominates. Moreover, the BFSS potential also has 
flat directions. Both of these observations combined could conceivably produce anomalous 
powers of N in the final answer, so it is worth checking that factorization holds. 

We will proceed in two steps. First, we will check some consequences of factorization at 



some large value of N. For example, consider the matrix model correlators (as in Eq. (23)) 
of the following form 

(O n L 0™) = (24) 

where Ol = tr(Z L ), and Ol = tr(Z L ), for the case L > 2 and we take Z = X 1 +iX 2 , or any 
of its rotations. Rotational invariance of the ensemble implies that A m ^ n = A m 5 m ^ n . Our 
goal is to understand the at large N. Notice that (Ol) = 0, so this is exactly one of 
the cases where the naive large iV factorization does not apply. Instead, we can consider A\ 
as our first non-trivial value, and use it to normalize the answers. We expect that because 
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the potential in the BFSS matrix model is a scaling function, that the ratios A^j A\ are 
independent of the effective coupling constant (3. Arguing analogously to [6], we can think 
of tr(Z L ) as a raising operator for a composite field (an in single string state), and tr(Z L ) as 
the corresponding lowering operator (an out single string state). The effective propagator for 
raising and lowering would be just A\. This is the naive argumentation if planar diagrams 
were applicable. Then we would find that 



A, 



m\(A 1 



(25) 



from all the free contractions between raising and lowering operators. This would be the 
leading diagram for closed string propagation without interactions, and furthermore, other 
diagrams with interactions would be suppressed by 1/N 2 . Thus the statistical distribution 
of tr(Z L ) would be that of a random Gaussian variable. A simple check is to see if this is 
correct is to bin the results of sampling the real part of tr(Z L ), divided by its normalization 
(which in this case is A^/y/2) and to compare it to a Gaussian model for the distribution 
normalized to the number of samples. This is independent of L. We can see the results of 



this procedure in Fig. 14 



Counts 




FIG. 14. Test of factorization for N = 87. We bin the samples of Re(C?£,(i)) for various L obtained 
from the time series after thermalization. We compare to a Gaussian model of the data. Using 
logarithmic scaling in the counts permits us to check the tails of the distribution. We put measured 
counts of zero at 10~ 4 . 

The results of the test are that the different Ol have Gaussian statistics and we conlcude 
that the correlators do factorize in this sense. We did this for iV = 87, but the results are 
very similar for other values of N. 

We also consider correlators such as 

( o L o M o L+M ) ^ gW+M + 0(1/jv3)] (26) 

^A\AfA L + M N 

which should give rise to structure constants C that have a well defined large N limit. If 
we ignore the 1/iV 3 corrections, then the C should be independent of N up to statistical 
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N = 10 


N = 13 


N = 18 


AT = 87 


C2,2,4 
63,3,6 


4.97 ±0.51 
6.97 ±0.86 


4.54 ±0.15 
6.9 ±0.4 


4.94 ±0.8 
7.58 ±0.5 


4.97 ± 1.2 
8.36 ± 1.4 



TABLE III. Values of Cl,m,l+m at various values of iV 



uncertainties. Rotation invariance implies that the correlators above are real, since Z can 
be obtained by an SO (9) rotation of Z. 

Note that this correlator decays only 1/N and not as 1/N 2 . This is important for un- 
derstanding the large error bars of the measurement The value of an instantaneous 
measurement on the left is of order one while the expectation value is of order 1/N, thus the 
various measurements must cancel each other most of the time, leaving a small residual. A 
non-zero average could be also called a "violation of Gaussianity" if we think of the Ol as 



statistically independent variables. A simple test for two possible C is shown in table III 



where we can see that the C are indeed A" independent given the error bars. This gives us 
confidence that the standard large A" counting is applicable. 

We can generalize Eq. (J26| to include time dependence and check that 



(0 L (t)0 M (t)0 L+M (t + a)) t C L 



L AM aL+M 



A{AfA{ 



N 



+ 0(1/N S 



(27) 



where now the Cl,m,l+m indicate nonlinear correlations with time dependence. If the matrix 
model and gravity are to be matched, these powers of A" should be robust. We expect that 
the right hand side will decay with time a, as correlations typically do in chaotic systems. 
A plot of the correlation function C2 } 2,&{ a )N~ l can be seen in Fig. 15, where it decays as 
expected. The statistical error band that one should associate to the graph is similar in size 
to the error bars seen in Table 




FIG. 15. Correlation function 6*2,2, &{a)/N compared to the normalized autocorrelations A 2 ' 4 (a) = 
(tr(Z 2 - 4 (t))tr(Z 2 ' 4 (t ± a))) t . The normalization is A 2 / 4 = A 2 { A {0) for L = 2,4 respectively. The 
statistical error band on 6*2,2,4 should be roughly at 25% of the maximum around the value at zero. 
The graph shows data from a run at N = 87. 

4 Similar issues appeared in [47], where an object similar to Cl,m,l+m had a theory prediction that was 
being tested against a model. 
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Note how the correlation between the different variables Oi{t) peaks when the auto- 
correlation of the variables also peaks and that the correlations between different variables 
roughly decay as the autocorrelation function decays. This is consistent with naive expec- 
tations. The purpose of this check is to show that large N counting is also applicable to 
general dynamical questions: if the initial non-gaussianity is of order 1/N and it bounds the 
time dependent non-gaussianity, then these can not be larger than 1/N. 



VII. GRAVITATIONAL INTERPRETATION 

In this section we discuss the extent to which we can call the thermalized classical dynam- 
ics of the matrix model a black hole. Let us begin with the D-brane background geometries 
in the absence of excitations. They are characterized by the supergravity solutions found in 
|48j . In the string frame the ten dimensional metric is given by 

ds 2 = #- 1/2 (r)(cfejj) + H l ' 2 {r)(dr 2 + r 2 dn 2 8 _ p ) } (28) 

where dx\\ are the p + 1 coordinates that run along the worldvolume of the p-brane, r is a 
radial direction and Qg-p is an 8 — p dimensional sphere, and if is a harmonic form given 
by 

H{r) = 1 + ^ (29) 

where a is a constant that depends on p, but not on N or r. 

We start with a few comments on string theory in this geometry. For p < 3 these metrics 
produce a long throat near r = 0, in that if 1//4 (r) dr = oo. If not for the time warping, 
a string that stretches to the brane, at the origin, would have infinite mass. This means 
the region near r = can be considered a large volume, in string units, and we are justified 



in dropping the 1 in (29). As noted in [7], the string coupling remains finite as we take 
the near horizon, decoupling limit. For p = 0, which is the case of interest given the SO(9) 
symmetry of the BFSS system, the effective curvature becomes large for large r, and is small 
near r = 0. 

Now imagine adding energy to this system and creating a black hole with the same asymp- 
totics as this background. These black holes have positive specific heat. If the temperature 
is low, the curvature of the black hole near the horizon is small in string units [7j . Because 
the dilaton runs in these geometries towards small coupling in the UV, the effective string 
scale depends on the position, and the curvature becomes large in string units in the UV. 
Correspondingly, the curvature will be large in string units if the temperature is high, so the 
high temperature black hole is stringy. The regime of interest for us is high temperature, 
where the curvature near the putative back hole is large and where stringy corrections must 
be important. In the spirit of [10], we should be able to describe that region by replacing the 
black holes by configurations of D-branes, because we are in the stringy regime. The classical 
dynamics we studied is the microscopic, classical dynamics of the D-brane configurations. 
These configurations with D-brane sources are horizonless in the sense of classical gravity, 
but when the system cools enough we recover a black hole. This notion that we recover a 
black hole when we cool the system down is the reverse of the scenario in [10J, where the 
authors argued that a black hole that becomes small enough to become stringy is replaced 
by a set of strings and possibly D-branes. We assume that this philosophy is applicable in 
this case as well, as there is no reason to expect a phase transition and we can corroborate 
this by examining the scalings of various quantities with respect to N and T. 
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From our calculations via the virial theorem, we know that the radius of the brane 
configurations in the classical system grows as R ~ jV 1//4 T 1//4 , so the effective (codimension 
two) size covered by the brane system scales as R s ~ N 2 T 2 . The energy is of order N 2 T 
whereas the entropy of the D-brane gas is of order iV 2 log (T/T ), where T is some reference 
temperature whose precise value doesn't concern us here. As computed in [33], for low 
temperatures the entropy scales as S ~ jV 2 T 9//5 and the free energy scales as£~ iV 2 T 14 / 5 . 
The smooth transition to the brane gas occurs near T ~ 1 in their units. Since for the 
black brane we have the usual area law for the entropy S ~ r 8 , for T ~ 1 we have R 8 ~ r 8 . 
So the brane gas extends all the way to where we would imagine the black hole horizon 
to be. The entropy in both cases is of order N 2 , as well as their energy, so there is no 
gap between the energy or entropy scalings that would suggest a phase transition. This is 
the essence of the argument in [TU] for smooth interpolation between strings and black hole 
physics. Furthermore, that the configuration reaches the same radius as where the black hole 
horizon would be located is very similar to the fuzzball geometries (see [50] for a review). 
This suggests that the system can be described both as a black hole and as a collection of 
D-branes. If the transition between the descriptions is smooth, we can choose one or the 
other depending on the question we want to ask. 

Now we want to consider the dissipation we observe, in the sense of decaying correlations, 



in our simulations. As shown in our Fig. 10, for all N we get a very similar power spectrum 
of fluctuations. Applying fluctuation-dissipation relations, we can use these to characterize 
response functions. If we were to match to gravity, valid at lower temperatures, we would 
expect dissipation to be related to the presence of quasinormal modes associated with a 
black hole horizon. 

We expect the point of comparison between both descriptions to be the analytic structure 



of correlation functions and power spectra, such as those shown in Fig. 10 As shown in 
the figure, there are many suggestive straight lines that describe the logarithm of the power 
spectrum. Remember also that the power spectrum is symmetric about u = 0. The most 
naive fit to the graph near u ~ is 

P(u) ~ exp(-/3|w|). (30) 
The absolute value is not an analytic function of the complex variable u and the Fourier 



transform of Eq. (30 ) only decays polynomially at large times. If the singularity is smoothed 
out very near u = 0, one can imagine that an analytic function of u replaces \u\. The simplest 
such function is f(uj) = Vu 2 ~+~e. When e — > we recover the result above. The presence 
of e suggests a pair of branch cuts beginning near u ~ along the imaginary axis. If these 
are the closest singularities to the real axis, then the function decays exponentially in time. 
Thus, even if the ultimate late time fate of the system has correlation functions that decay 
exponentially in time, there can be a long transient where the decays of the correlation 
functions are only polynomial. 

A square root branch cut can be approximated by a density of poles (this is often seen 
in matrix models [51]). Since we only have the analytic function on the real numbers, 
extrapolating to find the poles requires a very good understanding of the analytic structure 
of the function and the pattern of pole locations. Just for comparison, another similar 
function would be given by 

P{oo) — exp (-y/PuP + e- v //3*a; 2 + e*) , (31) 
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where /3, f3* are complex conjugates of each other. This is real on the real axis and has four 
branch cuts near the origin starting at u = ±iy/e//3 and u = ±zv/e*//3*. Such branch cuts 
could indicate a series of poles along straight lines, starting from the brach cut endpoints 
and going towards infinity. Such patterns of quasinormal modes aligning on fairly straight 
curves have been observed in Schwarzschild and AdS black holes A particularly 

nice set of examples can be found in [55], Fig. 4, and [56]. In gravity setups such patterns 
are interpreted in terms of the membrane paradigm [57] . 

Our results suggest that there is agreement between our analytics and the existence of 
a large collection of quasinormal modes for each L near u> = 0. We might eventually be 
able to interpret them as shear or sound modes once we understand the details of how these 
modes are mapped into each other. Notice that as we increase L, as in Fig. [12J the curves 
become flatter near u> = for even L, suggesting that the corresponding poles nearest to 
the origin are moving away from the real axis. This suggests that there is a dispersion 
relation for the frequencies uj(L) of these poles such that lm(u(L)) increases with L. Here 
L is angular momentum about the sphere of the spherical black hole geometry, so such a 
dispersion relation could be interpreted geometrically, but we haven't made this precise yet: 
in the absence of a theory of where these poles and branch cuts should be located, any match 
to a dispersion relation would be mere speculation. 

The analytic structure of the power spectra of observables in a chaotic dynamical system 
is generally controlled by poles off the real axis known as Pollicot-Ruelle resonances. The 
locations of these are known for some simple systems, but finding them in general, even 
for very low dimensional systems, is difficult (see, e.g., [58]). For large systems with many 
degrees of freedom there are expected to be accumulations of poles and also branch cuts 
[5*U] , in line with our discussion above. More detailed information on the analytic structure 
may be difficult to obtain, but we plan to analyze this in more detail in the future. We also 



note the other lines drawn in the Fig. 10, which suggest similar interpretations for the other 



peaks as modes whose frequencies begin with a non-zero real part. 

We do not have a theory for the analytic structure of the power spectra, but based on 
the data and our general understanding of dissipative phenomena we have speculated that 
there are branch cuts in the analytically continued spectra. Given the behavior of black hole 
horizons, as in the membrane paradigm and as characterized by the quasinormal modes, we 
think our numerical data are consistent with having a smooth transition from the matrix 
configurations to a black hole at low temperatures. 



VIII. CONCLUSION 

Let us first summarize our observations from our numerical simulations of the BMN 
and BFSS matrix models. The dynamics of the matrix models are chaotic. There are 
strong indications that the systems thermalize, most prominently that the typical matrices 
of the momentum and position variables behave, at late times, as random matrices from 
the traceless GUE ensemble. This is expected for the momentum matrices because the 
Hamiltonian is quadratic in the relevant degrees of freedom, but because of the nonlinearities 
in the potential we would not immediately guess this for the position matrices. We also 
showed how the presence of constraints alters the naive arguments about the appropriate 
random matrix ensembles for these systems. We have seen that certain observables behave 
hydrodynamically, i.e., their power spectra at equilibrium are approximately independent 
of the total number of degrees of freedom and they have approximately gaussian statistics. 
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We have also seen that large N counting applies to time correlation functions between 
observables, so that violations of gaussianity scale like 1/N, i.e., there is factorization of 
these degrees of freedom. This 1/N scaling is associated with quantum corrections under 
the usual AdS/CFT power counting arguments. 

These particular matrix models are interesting because of their dual gravitational inter- 
pretation. Though quantum effects are crucial to the emergent gravitational dynamics in 
general, we have argued that at high temperatures the classical dynamics of these matrix 
models do encode certain dual gravitational phenomena, including black holes, albeit in a 
stringy regime. As such, a direct comparison with supergravity solutions is not possible, but 
we can look for and have found some qualitative agreement. 

The first level of correspondence is between black holes, broadly defined to include large 
stringy corrections, and equilibrium thermal states in these models. This thermodynamic 
correspondence was referred to in both of the seminal papers on BFSS matrix theory and 
AdS/CFT [5], [22] and subsequently investigated and verified in many situations. The main 
work of this paper has been the identification of high temperature equilibrium configurations 
in the matrix models and some investigation of their equations of state and fluctuations, 
in the manner of a molecular dynamics simulation. This is a continuation of the work first 
presented in [13]. As a first step toward a more detailed correspondence, we have presented 
evidence here that the approximate analytic structure we find in the power spectra of certain 
observables at large N could be the remnant of the quasinormal modes of supergravity black 
holes. This is suggested by the numerical data, where the results we find on the real axis 
seem to be well approximated by functions that, when analytically continued, would have 
branch cuts in the complex plane. These could be an approximation of a sequence of roughly 
evenly spaced poles. 

It is possible to study the geometry that emerges from these matrix models in more 
detail, for example by using the methods presented recently in |60j. Such geometric data 
is based on probe branes and the fermionic degrees of freedom that connect them to the 
configurations we have shown here. This information would complement what is done here 
very nicely and we can speculate that in the future such reasoning will lead us to understand 
the emergence of gravitational horizons from such dynamics better. 

We also seek to investigate the dissipation and other transport properties further. For 
example, the time autocorrelation functions of various observables, such as those displayed 



in Fig. 15, can be integrated to compute the associated transport coefficients, via the 
Green-Kubo relations. We can then investigate how these depend on the temperature, N, 
or other variables, and see whether the behavior can be interpreted holographically and 
meets our expectations from gravity. We note that these transport phenomena are related 
to the chaotic dynamics of the system, characterized by quantities such as the Lyapunov 
exponents, Pollicot-Ruelle resonances, and Kolmogorov- Sinai entropy (see, e.g., [HI] for a 
review). These quantities also control the far from equilibrium dynamics. In that vein, 
we may be able to profitably study other non-equilibrium phenomena by applying further 
techniques and results from non-equilibrium statistical mechanics. 

Reflecting on what we have accomplished here, we have made a modest attempt at rec- 
onciling gravitational dynamics in holographic setups with the theory of chaotic dynamical 
systems. We have argued that these types of analyses do cover interesting black holes in the 
stringy regime, which would mean the study of these situations is forced on us if we want to 
be complete. From this vantage point, we have barely scratched the surface of what can be 
computed and analyzed and it would be interesting to pursue this further to improve our 
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understanding of holography. These simulations of the classical dynamics provide new ways 
to address questions about black holes in simple systems where the numerical computations 
are easily implemented. It is also important to understand how to move away from the 
correspondence limit, which would mean quantum dynamics starts becoming important and 
the dynamics of fermions starts affecting the results. In such setups the tools of quantum 
chaos will play an increasingly important role. This is closely related to the question of 
whether at low temperatures it is the fermionic degrees of freedom that dominate, or if they 
only play a marginal role in most of the dynamical regimes of interest. 
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Appendix A: Measuring the Temperature 

Here we derive the relation between the temperature of a system with constraints and the 
second moments of its degrees of freedom. Consider a system with 2D degrees of freedom 
x, p with the standard kinetic energy, potential V(x), and k constraints C l (x,p) = with 
each C l linear in the momenta. The canonical partition function is given by 



Z = f d D xd D p Y[5(Ci(x,p))exp 

i=i 



(Al) 



We scale the momenta by a parameter yfiy and then rescale the 5 functions by the inverse 
of that parameter. Since the constraints are linear in the momenta 



i=l 



Z = 7 ( D - fe )/ 2 / d D xd D p U<f(Ci(x,p))exp -13 (J-\p\ 2 + V(x 



(A2) 
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The partition function is independent of 7 since all we have done is rescaled the momenta. 
Differentiating we have 



n = — = D ~ k ~,(D-k)/2 



d D xd D p l[5(C t (x,p))exp [-(3 (||pf + V(x)J 

f rf D ^ D pb1 2 n5(a(f, P 1)exp[-/3(||pl 2 + \/(f))] (A3) 
■' i=i 



Letting 7 = 1 we have 



= °-^Z - P -Z(\p\ 2 ) ^{D- k)T = (|pf ) 



(A4) 



In order to measure the temperature properly, we must subtract off the number of constraints 
from the degrees of freedom. 
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